###########################################################################
# IMPORTANT: Run via '00_master.R' or ensure WD is 'Codes' folder
###########################################################################

# --- 1. Path Management (Replication Safe) ---
# If this script is run standalone (not from Master), define defaults.
# This ensures it works if you just open this file and click "Run".
if (!exists("tables_dir")) tables_dir <- "Tables"
if (!exists("data_dir"))   data_dir   <- "Data"


# Create Tables folder if running standalone and it's missing
if (!dir.exists(tables_dir)) dir.create(tables_dir, showWarnings = FALSE)

# --- 2. Load Data ---
# Full sample
analysis <- read_dta(file.path(data_dir, "analysis.dta"))

# Factor labeling
analysis$conditions <- factor(analysis$conditions, labels = c("Control", "Placebo", "Treatment", "Existing"))  
analysis$spectrum <- factor(analysis$spectrum, labels = c("Pan-Blue", "Nonpartisan", "Pan-Green")) 

# Non-experimental sample
existing.df <- analysis %>%
  filter(conditions== "Existing") 

# Experimental sample
experiment.df <- analysis %>%
  filter(conditions!= "Existing") 


###################################
############ Main Text ############ 
###################################

# Table 1: Sample: Summary statistics, attrition, and balance test
demographics <- c("age", "female", "special_municipality", "edu", "full_time", 
                  "married", "income_cat", "pan_blue", "independent", "pan_green", 
                  "president_2016_voted", "president_2012_voted", "father_islander", 
                  "taoism", "nwpaper_pref_apple", "sanlih_show") #demographic controls
all_vars <- c(demographics, "control", "placebo", "treatment")
results <- list()

for (X in all_vars) {
  # Column 1: baseline survey among all subjects
  col1 <- analysis[[X]]
  N_w1 <- sum(!is.na(col1))
  mu_w1 <- mean(col1, na.rm = TRUE)
  sd_w1 <- sd(col1, na.rm = TRUE)
  # Column 2: endline survey among all subjects
  col2 <- analysis[[X]][analysis$wave2 == 1]
  N_w2 <- sum(!is.na(col2))
  mu_w2 <- mean(col2, na.rm = TRUE)
  sd_w2 <- sd(col2, na.rm = TRUE)
  # Column 3: attrition test (diff-attrition-of-exp-subjects)
  subset_attrition <- analysis[analysis$china_times != 1, ]
  ttest_result <- t.test(subset_attrition[[X]] ~ subset_attrition$wave2, var.equal = TRUE)
  attrition_pvalue <- ttest_result$p.value
  # Column 4: CT group in endline
  col4 <- analysis[[X]][analysis$china_times == 1 & analysis$wave2 == 1]
  N_ct <- sum(!is.na(col4))
  mu_ct <- mean(col4, na.rm = TRUE)
  sd_ct <- sd(col4, na.rm = TRUE)
  # Column 5: Control group in endline
  col5 <- analysis[[X]][analysis$control == 1 & analysis$wave2 == 1]
  N_c <- sum(!is.na(col5))
  mu_c <- mean(col5, na.rm = TRUE)
  sd_c <- sd(col5, na.rm = TRUE)
  # Column 6: Placebo group in endline
  col6 <- analysis[[X]][analysis$placebo == 1 & analysis$wave2 == 1]
  N_p <- sum(!is.na(col6))
  mu_p <- mean(col6, na.rm = TRUE)
  sd_p <- sd(col6, na.rm = TRUE)
  # Column 7: Treatment group in endline
  col7 <- analysis[[X]][analysis$treatment == 1 & analysis$wave2 == 1]
  N_t <- sum(!is.na(col7))
  mu_t <- mean(col7, na.rm = TRUE)
  sd_t <- sd(col7, na.rm = TRUE)
  # Column 8: ANOVA of experimental subjects
  subset_anova <- analysis[analysis$wave2 == 1 & analysis$china_times != 1, ]
  anova_result <- aov(as.formula(paste(X, "~ as.factor(conditions)")), data = subset_anova)
  anova_summary <- summary(anova_result)
  anova_pvalue <- anova_summary[[1]][["Pr(>F)"]][1]
  # Store results
  results[[X]] <- data.frame(
    variable = X,
    mu_w1 = mu_w1, sd_w1 = sd_w1, N_w1 = N_w1,
    mu_w2 = mu_w2, sd_w2 = sd_w2, N_w2 = N_w2,
    attrition_pvalue = attrition_pvalue,
    mu_ct = mu_ct, sd_ct = sd_ct, N_ct = N_ct,
    mu_c = mu_c, sd_c = sd_c, N_c = N_c,
    mu_p = mu_p, sd_p = sd_p, N_p = N_p,
    mu_t = mu_t, sd_t = sd_t, N_t = N_t,
    anova_pvalue = anova_pvalue)}

summstats_main <- bind_rows(results)
format_cols <- c("mu_w1", "sd_w1", "mu_w2", "sd_w2", "attrition_pvalue",
                 "mu_ct", "sd_ct", "mu_c", "sd_c", "mu_p", "sd_p", 
                 "mu_t", "sd_t", "anova_pvalue")
summstats_main[format_cols] <- lapply(summstats_main[format_cols], 
                                      function(x) sprintf("%.3f", x))

stargazer(summstats_main, 
          type = "latex",
          out = file.path(tables_dir, "Table 1.tex"), # Updated Output Path
          title = "Table 1: Attrition and Balance Test",
          summary = FALSE, rownames = FALSE, header = FALSE)


# Table 2: Browsing behavior and compliance

# Panel A: For all participants

#Prepare the Data
browsing_summary_df <- experiment.df %>%
  filter(treatment == 1 | placebo == 1) %>%
  mutate(condition = if_else(treatment == 1, "Treatment", "Placebo"),
         is_visited = if_else(!is.na(avg_browsing), 1, 0),
         is_min      = if_else(!is.na(avg_browsing) & avg_browsing >= 60, 1, 0),
         is_full     = if_else(!is.na(avg_browsing) & avg_browsing >= 180, 1, 0)
  ) %>%
  group_by(condition) %>%
  summarise(
    # Visited
    `Visited (Prop)` = mean(is_visited),
    `Visited (SD)`   = sd(is_visited),
    
    # Minimum >= 60s
    `>= 60s (Prop)`  = mean(is_min),
    `>= 60s (SD)`    = sd(is_min),
    
    # Full >= 180s
    `>= 180s (Prop)` = mean(is_full),
    `>= 180s (SD)`   = sd(is_full)
  )

#Reshape for the Table production
browsing_tableA <- browsing_summary_df %>%
  pivot_longer(cols = -condition, names_to = "Metric", values_to = "Value") %>%
  mutate(Value = sprintf("%.4f", Value)) %>%
  pivot_wider(names_from = condition, values_from = Value) %>%
  # Set column order
  select(Metric, Treatment, Placebo) %>% 
  # Set row order
  mutate(Metric = factor(Metric, levels = c(
    "Visited (Prop)", "Visited (SD)",
    ">= 60s (Prop)", ">= 60s (SD)",
    ">= 180s (Prop)", ">= 180s (SD)"))) %>%
  arrange(Metric) %>%
  mutate(Metric = as.character(Metric))


#Output
stargazer(browsing_tableA, 
          type = "latex", 
          out = file.path(tables_dir, "Table 2A.tex"),
          title = "Panel A: For all participants",
          summary = FALSE, rownames = FALSE, header = FALSE)


#Panel B
browsing_summary_df2 <- experiment.df %>%
  filter(treatment == 1 | placebo == 1) %>%
  mutate(condition = if_else(treatment == 1, "Treatment", "Placebo"),
         browsing_mins = avg_browsing / 60) %>%
  group_by(condition) %>%
  summarise(
    #Visited
    `Visited (Mean)` = mean(browsing_mins, na.rm = TRUE),
    `Visited (SD)`   = sd(browsing_mins, na.rm = TRUE),
    # Minimum >= 60s
    `>= 60s (Mean)`  = mean(browsing_mins[browsing_mins >= 1], na.rm = TRUE),
    `>= 60s (SD)`    = sd(browsing_mins[browsing_mins >= 1], na.rm = TRUE),
    # Full >= 180 seconds
    `>= 180s (Mean)` = mean(browsing_mins[browsing_mins >= 3], na.rm = TRUE),
    `>= 180s (SD)`   = sd(browsing_mins[browsing_mins >= 3], na.rm = TRUE))


#Reshape
browsing_tableB <- browsing_summary_df2 %>%
  pivot_longer(cols = -condition, names_to = "Metric", values_to = "Value") %>%
  mutate(Value = sprintf("%.2f", Value)) %>% 
  pivot_wider(names_from = condition, values_from = Value) %>%
  select(Metric, Treatment, Placebo) %>%     
  mutate(Metric = factor(Metric, levels = c(
    "Visited (Mean)", "Visited (SD)",
    ">= 60s (Mean)", ">= 60s (SD)",
    ">= 180s (Mean)", ">= 180s (SD)"
  ))) %>%
  arrange(Metric) %>%
  mutate(Metric = as.character(Metric))

#Output
stargazer(browsing_tableB, 
          type = "latex", 
          out = file.path(tables_dir, "Table 2B.tex"),
          title = "Panel B: Average browsing time",
          summary = FALSE, rownames = FALSE, header = FALSE)



########################################################
############ Online Supplementary Materials ############ 
########################################################

## Table SI-1: Regression of PRC index on personal characteristics
table_si1 <- lm_robust(prc_index_w1 ~ pan_blue + pan_green + age + female + 
                         special_municipality + edu + full_time + married + 
                         income_cat + president_2016_voted + president_2012_voted + 
                         father_islander + taoism + nwpaper_pref_apple + sanlih_show, 
                       data = analysis)
#output
texreg(table_si1,
       file = file.path(tables_dir, "Table SI-1.tex"), # Updated Output Path
       title = "Table SI-1: Baseline PRC Index Scores and Demographics",
       omit.coef = "Intercept",
       include.ci = FALSE,
       digits = 2,
       custom.coef.names = c("Pan Blue", "Pan Green", "Age", "Female", 
                             "Residence", "Education", "Full-Time Job", "Married", 
                             "Income", "Vote in 2016", "Vote in 2012", 
                             "Ethnicity", "Taoism", "Apple Daily", "SET News Program"))



## Table SI-2: Outcome variables descriptive statistics
key.outcomes <- c("vote_han_w1", "vote_han_w2", "vote_han_change", 
                  "relative_thermo_w1", "relative_thermo_w2", "relative_thermo_change", 
                  "prc_index_w1", "prc_index_w2", "prc_index_change", 
                  "pan_blue", "independent", "pan_green")
var_labels <- c("Vote for Han (Baseline)", "Vote for Han (Endline)", "Vote for Han (Change Score)",
                "Candidate Evaluation (Baseline)", "Candidate Evaluation (Endline)", "Candidate Evaluation (Change Score)", 
                "Pro-PRC Index (Baseline)", "Pro-PRC Index (Endline)", "Pro-PRC Index (Change Score)",
                "PRC-friendly", "Nonpartisan", "PRC-skeptical")
all_obs <- as.data.frame(analysis[, key.outcomes])
exp_obs <- as.data.frame(experiment.df[, key.outcomes])
nonexp_obs <- as.data.frame(existing.df[, key.outcomes])

#output Panel A: All participants
stargazer(all_obs,
          type = "latex",
          out = file.path(tables_dir, "Table SI-2A.tex"), # Updated Output Path
          summary.stat = c("n", "mean", "sd", "min", "max"),
          digits = 2,
          covariate.labels = var_labels,
          title = "Table SI-2: All Participants",
          header = FALSE)
#output Panel B: Experimental participants
stargazer(exp_obs,
          type = "latex",
          out = file.path(tables_dir, "Table SI-2B.tex"), # Updated Output Path
          summary.stat = c("n", "mean", "sd", "min", "max"),
          digits = 2,
          covariate.labels = var_labels,
          title = "Table SI-2: Experimental Participants",
          header = FALSE)

#output Panel C: Non-experimental participants
stargazer(nonexp_obs,
          type = "latex",
          out = file.path(tables_dir, "Table SI-2C.tex"), # Updated Output Path
          summary.stat = c("n", "mean", "sd", "min", "max"),
          digits = 2,
          covariate.labels = var_labels,
          title = "Table SI-2: Nonexperimental Participants",
          header = FALSE)


## Table SI-3: Treatment effects on outcomes
vote.han.existing <- lm_robust(vote_han_change ~ 1, data = existing.df)
vote.han.itt <- lm_robust(vote_han_change ~ treatment, data = experiment.df)
vote.han.min <- iv_robust(vote_han_change ~ min_complier | treatment, data = experiment.df)
vote.han.full <- iv_robust(vote_han_change ~ full_complier | treatment, data = experiment.df)

favhan.existing <- lm_robust(relative_thermo_change ~ 1, data = existing.df)
favhan.itt <- lm_robust(relative_thermo_change ~ treatment, data = experiment.df)
favhan.min <- iv_robust(relative_thermo_change ~ min_complier | treatment, data = experiment.df)
favhan.full <- iv_robust(relative_thermo_change ~ full_complier | treatment, data = experiment.df)

china.existing <- lm_robust(prc_index_change ~ 1, data = existing.df)
china.itt <- lm_robust(prc_index_change ~ treatment, data = experiment.df)
china.min <- iv_robust(prc_index_change ~ min_complier | treatment, data = experiment.df)
china.full <- iv_robust(prc_index_change ~ full_complier | treatment, data = experiment.df)

#output Panel A and Panel B
texreg(list(vote.han.itt, vote.han.min, vote.han.full,
            favhan.itt, favhan.min, favhan.full,
            china.itt, china.min, china.full),
       include.ci = FALSE, include.rsquared = FALSE, include.adjrs = FALSE, include.rmse = FALSE,
       file = file.path(tables_dir, "Table SI-3AB.tex"), # Updated Output Path
       title = "Table SI-3: Panel A and B",
       digits = 3,
       stars = c(0.01, 0.05, 0.1),
       omit.coef = "Intercept",
       custom.coef.names = c("ITT", 
                             "TOT (Minimum)", 
                             "TOT (Full)"),
       custom.model.names = c("(1)", "(2)", "(3)", 
                              "(4)", "(5)", "(6)", 
                              "(7)", "(8)", "(9)"),
       custom.header = list("Vote for Han" = 1:3, 
                            "Candidate Evaluation" = 4:6, 
                            "Pro-PRC Index" = 7:9))

pooled_exp <- analysis %>%
  filter(conditions != "Existing" & wave2==1) %>%
  summarise(
    conditions = "Experimental Sample",
    mean_vote = mean(vote_han_change, na.rm = TRUE),
    sd_vote = sd(vote_han_change, na.rm = TRUE),
    mean_thermo = mean(relative_thermo_change, na.rm = TRUE),
    sd_thermo = sd(relative_thermo_change, na.rm = TRUE),
    mean_prc = mean(prc_index_change, na.rm = TRUE),
    sd_prc = sd(prc_index_change, na.rm = TRUE))

by_conditions <- analysis %>%
  filter(wave2==1) %>%
  group_by(conditions) %>%
  summarise(
    mean_vote = mean(vote_han_change, na.rm = TRUE),
    sd_vote = sd(vote_han_change, na.rm = TRUE),
    mean_thermo = mean(relative_thermo_change, na.rm = TRUE),
    sd_thermo = sd(relative_thermo_change, na.rm = TRUE),
    mean_prc = mean(prc_index_change, na.rm = TRUE),
    sd_prc = sd(prc_index_change, na.rm = TRUE))

summary_stats <- bind_rows(pooled_exp, by_conditions)

#output Panel C
stargazer(as.data.frame(summary_stats),
          type = "latex",
          out = file.path(tables_dir, "Table SI-3C.tex"), # Updated Output Path
          title = "Table SI-3: Panel C",
          digits = 3,
          summary = FALSE,
          header = FALSE)


## Table SI-4: Treatment effects on vote choice by political attentiveness
vote_engagement <- interflex(Y = "vote_han_change", D = "treatment", 
                             X = "concern_election_outcome_w1", data = as.data.frame(experiment.df), 
                             na.rm = TRUE, estimator = "linear", vcov.type = "robust", Xdistr = "density")

vote_engagement.estimate <- data.frame(vote_engagement$est.lin) %>%
  rename(
    "Moderator" = X1.X,
    "Marginal Effect" = X1.ME,
    "Robust SE" = X1.sd,
    "Lower CI" = X1.lower.CI.95..,
    "Upper CI" = X1.upper.CI.95..)

stargazer(as.data.frame(vote_engagement.estimate),
          type = "latex",
          out = file.path(tables_dir, "Table SI-4.tex"), # Updated Output Path
          title = "Table SI-4: Vote Choice",
          align = TRUE, summary = FALSE, rownames = FALSE, header = FALSE,
          digits = 3)

## Table SI-5: Treatment effects on candidate evaluation by political attentiveness
thermo_engagement <- interflex(Y = "relative_thermo_change", D = "treatment", 
                               X = "concern_election_outcome_w1", data = as.data.frame(experiment.df), 
                               na.rm = TRUE, estimator = "linear", vcov.type = "robust", Xdistr = "density")

thermo_engagement.estimate <- data.frame(thermo_engagement$est.lin) %>%
  rename(
    "Moderator" = X1.X,
    "Marginal Effect" = X1.ME,
    "Robust SE" = X1.sd,
    "Lower CI" = X1.lower.CI.95..,
    "Upper CI" = X1.upper.CI.95..)

stargazer(as.data.frame(thermo_engagement.estimate),
          type = "latex",
          out = file.path(tables_dir, "Table SI-5.tex"), # Updated Output Path
          title = "Table SI-5: Candidate Evaluation",
          align = TRUE, summary = FALSE, rownames = FALSE, header = FALSE,
          digits = 3)


## Table SI-6: Treatment effects on Pro-PRC index by political attentiveness
china_engagement <- interflex(Y = "prc_index_change", D = "treatment", 
                              X = "concern_election_outcome_w1", data = as.data.frame(experiment.df), 
                              na.rm = TRUE, estimator = "linear", vcov.type = "robust", Xdistr = "density")

china_engagement.estimate <- data.frame(china_engagement$est.lin) %>%
  rename(
    "Moderator" = X1.X,
    "Marginal Effect" = X1.ME,
    "Robust SE" = X1.sd,
    "Lower CI" = X1.lower.CI.95..,
    "Upper CI" = X1.upper.CI.95..)

stargazer(as.data.frame(china_engagement.estimate),
          type = "latex",
          out = file.path(tables_dir, "Table SI-6.tex"), # Updated Output Path
          title = "Table SI-6: Pro-PRC Index",
          align = TRUE, summary = FALSE, rownames = FALSE, header = FALSE,
          digits = 3)


## Table SI-7: Treatment effect by political predispositions
# PRC-friendly respondents
bluehan.itt <- lm_robust(vote_han_change ~ treatment, 
                         data = subset(experiment.df, spectrum=="Pan-Blue"))
bluehan.att <- iv_robust(vote_han_change ~ full_complier | treatment, 
                         data = subset(experiment.df, spectrum=="Pan-Blue"))
bluefavhan.itt <- lm_robust(relative_thermo_change ~ treatment, 
                            data = subset(experiment.df, spectrum=="Pan-Blue"))
bluefavhan.att <- iv_robust(relative_thermo_change ~ full_complier | treatment, 
                            data = subset(experiment.df, spectrum=="Pan-Blue"))
bluechina.itt <- lm_robust(prc_index_change ~ treatment, 
                           data = subset(experiment.df, spectrum=="Pan-Blue"))
bluechina.att <- iv_robust(prc_index_change ~ full_complier | treatment, 
                           data = subset(experiment.df, spectrum=="Pan-Blue"))
# Nonpartisan respondents
indephan.itt <- lm_robust(vote_han_change ~ treatment, 
                          data = subset(experiment.df, spectrum=="Nonpartisan"))
indephan.att <- iv_robust(vote_han_change ~ full_complier | treatment, 
                          data = subset(experiment.df, spectrum=="Nonpartisan"))
indepfavhan.itt <- lm_robust(relative_thermo_change ~ treatment, 
                             data = subset(experiment.df, spectrum=="Nonpartisan"))
indepfavhan.att <- iv_robust(relative_thermo_change ~ full_complier | treatment, 
                             data = subset(experiment.df, spectrum=="Nonpartisan"))
indepchina.itt <- lm_robust(prc_index_change ~ treatment, 
                            data = subset(experiment.df, spectrum=="Nonpartisan"))
indepchina.att <- iv_robust(prc_index_change ~ full_complier | treatment, 
                            data = subset(experiment.df, spectrum=="Nonpartisan"))
# PRC-skeptical respondents
greenhan.itt <- lm_robust(vote_han_change ~ treatment, 
                          data = subset(experiment.df, spectrum=="Pan-Green"))
greenhan.att <- iv_robust(vote_han_change ~ full_complier | treatment, 
                          data = subset(experiment.df, spectrum=="Pan-Green"))
greenfavhan.itt <- lm_robust(relative_thermo_change ~ treatment, 
                             data = subset(experiment.df, spectrum=="Pan-Green"))
greenfavhan.att <- iv_robust(relative_thermo_change ~ full_complier | treatment, 
                             data = subset(experiment.df, spectrum=="Pan-Green"))
greenchina.itt <- lm_robust(prc_index_change ~ treatment, 
                            data = subset(experiment.df, spectrum=="Pan-Green"))
greenchina.att <- iv_robust(prc_index_change ~ full_complier | treatment, 
                            data = subset(experiment.df, spectrum=="Pan-Green"))

# Updated Output Paths in list
subgroups <- list(
  list(name = "PRC-Friendly", file = file.path(tables_dir, "Table SI-7A.tex"),
       models = list(bluehan.itt, bluehan.att, bluefavhan.itt, bluefavhan.att, bluechina.itt, bluechina.att)),
  list(name = "Nonpartisan", file = file.path(tables_dir, "Table SI-7B.tex"),
       models = list(indephan.itt, indephan.att, indepfavhan.itt, indepfavhan.att, indepchina.itt, indepchina.att)),
  list(name = "PRC-Skeptical", file = file.path(tables_dir, "Table SI-7C.tex"),
       models = list(greenhan.itt, greenhan.att, greenfavhan.itt, greenfavhan.att, greenchina.itt, greenchina.att)))

make_subgroup_table <- function(subgroup) {
  texreg(subgroup$models,
         file = subgroup$file,
         digits = 3, omit.coef = "Intercept", label = "",
         include.ci = FALSE, include.rsquared = FALSE, include.adjrs = FALSE, include.rmse = FALSE,
         custom.model.names = c("ITT", "TOT", "ITT", "TOT", "ITT", "TOT"),
         custom.header = list("Vote for Han" = 1:2, "Candidate Eval" = 3:4, "Pro-PRC" = 5:6),
         caption = paste(subgroup$name, "Subgroup"))}

lapply(subgroups, make_subgroup_table)


## Table SI-8: Baseline vote intent among all experimental participants
crosstab <- analysis %>%
  filter(conditions != "Existing") %>%
  mutate(pres_vote_w1 = factor(pres_vote_w1, 
                               levels = c(1, 2, 3, 4),
                               labels = c("Han", "Tsai", "Soong", "Undecided"))) %>%
  count(spectrum, pres_vote_w1) %>%
  group_by(spectrum) %>%
  mutate(row_total = sum(n),
         cell = sprintf("%d (%.2f%%)", n, n / sum(n) * 100)) %>%
  select(-n) %>%
  pivot_wider(names_from = pres_vote_w1, values_from = cell, values_fill = "0 (0.0%)") %>%
  relocate(row_total, .after = last_col()) %>%
  rename(`Political Predisposition` = spectrum,
         `N` = row_total)

stargazer(as.data.frame(crosstab),
          type = "latex", out = file.path(tables_dir, "Table SI-8.tex"), # Updated Output Path
          summary = FALSE, rownames = FALSE, header = FALSE,
          title = "Vote Intent of Experimental Participants", label = "")


## Table SI-9: Baseline vote intent among experimental participants completing both waves
crosstab_both <- analysis %>%
  filter(conditions != "Existing" & wave2 == 1) %>%
  mutate(pres_vote_w1 = factor(pres_vote_w1, 
                               levels = c(1, 2, 3, 4),
                               labels = c("Han", "Tsai", "Soong", "Undecided"))) %>%
  count(spectrum, pres_vote_w1) %>%
  group_by(spectrum) %>%
  mutate(row_total = sum(n),
         cell = sprintf("%d (%.2f%%)", n, n / sum(n) * 100)) %>%
  select(-n) %>%
  pivot_wider(names_from = pres_vote_w1, values_from = cell, values_fill = "0 (0.0%)") %>%
  relocate(row_total, .after = last_col()) %>%
  rename(`Political Predisposition` = spectrum,
         `N` = row_total)

stargazer(as.data.frame(crosstab_both),
          type = "latex", out = file.path(tables_dir, "Table SI-9.tex"), # Updated Output Path
          summary = FALSE, rownames = FALSE, header = FALSE,
          title = "Vote Intent of Experimental Participants Completing Both Waves", label = "")


## Table SI-10: Regression of red-media perceptions on personal characteristics
red_media <- lm_robust(reformulate(demographics[demographics!="pan_green"], "ct_rmedia"), 
                       data = subset(analysis, spectrum!="Nonpartisan"))
texreg(red_media,
       file = file.path(tables_dir, "Table SI-10.tex"), # Updated Output Path
       digits = 3, omit.coef = "Intercept",
       include.ci = FALSE, include.rsquared = FALSE, include.adjrs = FALSE, include.rmse = FALSE, header = FALSE)


## Table SI-11: Regression of news consumption on treatment conditions
news_consumption <- lm_robust(reformulate(c("treatment", demographics[demographics!="independent"]), 
                                          "overall_polinews_consumption_w2"), data = analysis)
election_news <- lm_robust(reformulate(c("treatment", demographics[demographics!="independent"]), 
                                       "pnews_days_w2"), data = analysis)

texreg(list(news_consumption, election_news),
       file = file.path(tables_dir, "Table SI-11.tex"), # Updated Output Path
       digits = 3,
       omit.coef = "Intercept",
       custom.model.names = c("News Consumption", "Election News"),
       custom.coef.names = c("Treatment", "Age", "Female", "Residence", "Education",
                             "Full-Time Job", "Married", "Income", "Pan-Blue", 
                             "Pan-Green", "Vote in 2016", "Vote in 2012", 
                             "Ethnicity", "Taoism", "Apple Daily", "SET News Program"),
       include.ci = FALSE, include.rsquared = FALSE, include.adjrs = FALSE, include.rmse = FALSE, header = FALSE,
       caption = "Table SI-11: News Consumption and Treatment Conditions")


## Table SI-12: Covariate-adjusted estimates
vote_covariate <- lm_lin(vote_han_change ~ treatment, 
                         covariates = ~ pan_green + pan_blue + age + female + special_municipality +
                           edu + full_time + married + income_cat + president_2016_voted + president_2012_voted +
                           father_islander + taoism + nwpaper_pref_apple + sanlih_show, data = experiment.df)
favhan_covariate <- lm_lin(relative_thermo_change ~ treatment, 
                           covariates = ~ pan_green + pan_blue + age + female + special_municipality +
                             edu + full_time + married + income_cat + president_2016_voted + president_2012_voted +
                             father_islander + taoism + nwpaper_pref_apple + sanlih_show, data = experiment.df)
china_covariate <- lm_lin(prc_index_change ~ treatment, 
                          covariates = ~ pan_green + pan_blue + age + female + special_municipality +
                            edu + full_time + married + income_cat + president_2016_voted + president_2012_voted +
                            father_islander + taoism + nwpaper_pref_apple + sanlih_show, data = experiment.df)

texreg(list(vote_covariate, favhan_covariate, china_covariate),
       file = file.path(tables_dir, "Table SI-12.tex"), # Updated Output Path
       digits = 2,
       omit.coef = "Intercept",
       custom.model.names = c("Vote for Han", "Candidate Evaluation", "Pro-PRC Index"),
       custom.coef.names = c("Treatment",
                             "Pan-Green", "Pan-Blue", "Age", "Female", "Residence",
                             "Education", "Full-Time Job", "Married", "Income", 
                             "Vote in 2016", "Vote in 2012", "Ethnicity", "Taoism", 
                             "Apple Daily", "SET News Program",
                             "Treatment × Pan-Green", "Treatment × Pan-Blue", 
                             "Treatment × Age", "Treatment × Female", "Treatment × Residence",
                             "Treatment × Education", "Treatment × Full-Time Job", 
                             "Treatment × Married", "Treatment × Income", 
                             "Treatment × Vote in 2016", "Treatment × Vote in 2012", 
                             "Treatment × Ethnicity", "Treatment × Taoism", 
                             "Treatment × Apple Daily", "Treatment × SET News Program"),
       include.ci = FALSE, include.rsquared = FALSE, include.adjrs = FALSE, include.rmse = FALSE, header = FALSE,
       caption = "Table SI-12: Covariate-Adjusted ITT Estimates")


## Table SI-13: multiple comparisons and adjusted p-values 
p_values <- c(
  summary(vote.han.itt)$coefficients[2, 4],
  summary(bluehan.itt)$coefficients[2, 4],
  summary(indephan.itt)$coefficients[2, 4],
  summary(greenhan.itt)$coefficients[2, 4],
  summary(favhan.itt)$coefficients[2, 4],
  summary(bluefavhan.itt)$coefficients[2, 4],
  summary(indepfavhan.itt)$coefficients[2, 4],
  summary(greenfavhan.itt)$coefficients[2, 4],
  summary(china.itt)$coefficients[2, 4],
  summary(bluechina.itt)$coefficients[2, 4],
  summary(indepchina.itt)$coefficients[2, 4],
  summary(greenchina.itt)$coefficients[2, 4])


p_adjust_bonferroni <- p.adjust(p_values, method = "bonferroni") # Bonferroni
p_adjust_fdr <- p.adjust(p_values, method = "fdr") # FDR correction

p_adjusted_results <- data.frame(
  Outcome = c("Vote for Han", 
              "Vote for Han (PRC-friendly)", 
              "Vote for Han (Nonpartisan)", 
              "Vote for Han (PRC-skeptical)",
              "Candidate Eval.", 
              "Candidate Eval. (PRC-friendly)",
              "Candidate Eval. (Nonpartisan)", 
              "Candidate Eval. (PRC-skeptical)",
              "PRC Index", 
              "PRC Index (PRC-friendly)", 
              "PRC Index (Nonpartisan)", 
              "PRC Index (PRC-skeptical)"),
  P_Value = p_values,
  Bonferroni_Corrected = p_adjust_bonferroni,
  FDR_Corrected = p_adjust_fdr) %>%
  stargazer(type = "latex",
            out = file.path(tables_dir, "Table SI-13.tex"), # Updated Output Path
            title = "Table SI-13: Adjusted P-Values",
            digits = 6,
            summary = FALSE,
            header = FALSE)


cat("Tables successfully reproduced. End of table-replication script.\n")